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ABSTRACT 

The effects of initially uniform magnetic fields on the formation and evolution of dense pillars 
and cometary globules at the boundaries of H II regions are investigated using 3D radiation- 
magnetohydrodynamics simulations. It is shown, in agreement with previous work, that a 
strong initial magnetic field is required to significantly alter the non-magnetised dynamics 
because the energy input from photoionisation is so large that it remains the dominant driver 
of the dynamics in most situations. Additionally it is found that for weak and medium field 
strengths an initially perpendicular field is swept into alignment with the pillar during its dy- 
namical evolution, matching magnetic field observations of the 'Pillars of Creation' in Ml 6 
and also some cometary globules. A strong perpendicular magnetic field remains in its ini- 
tial configuration and also confines the photoevaporation flow into a bar- shaped dense ionised 
ribbon which partially shields the ionisation front and would be readily observable in recom- 
bination lines. A simple analytic model is presented to explain the properties of this bright 
linear structure. These results show that magnetic field strengths in star-forming regions can 
in principle be significantly constrained by the morphology of structures which form at the 
borders of H II regions. 
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1 INTRODUCTION 

The study of elephant trunks, pillars and globules, found at the 
borders of H II regions around massive stars, has received sig- 
nificant attention in recent years, both from an observational per- 
spective and in theoretical and computational models. The fa- 
mous 'Pillars of C reation' in M l 6 were observe d at optical wave- 
lengths with HST (iHester et aLl[T99 6). in the IR (Indebetouw et al. 
2007l:ISugitani et al.l2007h . and in sub-mm/radio (e.g. iPound 1998. : 
White etal.ll 19991) . showing that these are dynamic structures with 
ongoing star formation which may or may not have been triggered 
by the radiation which has shaped the pillars. Smith et al. ( 2010b) 
provide strong evidence that pillars in the Carina Nebula are signif- 
icant sites of sequential star formation propagating away from the 
older star clusters in this region, building on previous observations 
of synchronised star formation around the periphery of the nebula 
by Smith & Brooks (2007). On s maller scales , studies of T -Taur i 
star ages in the Orion Nebula Lee et^ [20051 : Lee & Chenll2007h 
show decreasing stellar ages moving away from massive stars and 
towards bright-rimmed clouds at the H II region/molecular cloud 
interface, again strongly suggesting at least sequential and possi- 
bly triggered star formation. The clear relationship between pil- 
lars/globules and second generation star formation around OB as- 
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sociations, and the question of the extent to which this star for- 
mation is triggered, provides strong motivation to understand the 
formation and evolution of these structures. 

In a previous paper (Mackev &Limll2010 '. hereafter MLIO) 
we investigated the formation and evolution of dense pillars of gas 
and dust - elephant trunks - on the boundaries of H II regions using 
3D hydrodynamical simulations including photoionising radiative 
transfer (R-HD). It was found that shadowing of ionising radia- 
tion by an inhomogeneous density field naturally forms elephant 
trunks without the assistance of self-gravity, or of ionisation front 

and cooling instabilities. A combination of radiation-driven implo- 
' in 1 

sion (RDI; Bertoldi 1989) and acceleration due to the rocket ef- 
fect ([OorL& Spitzer 1955) produce elongated structures: RDI com- 
presses neutral gas until pressure equilibrium with ionised gas is 
achieved; the rocket effect accelerates gas away from the radiation 
source producing dynamic elongated structures with lifetimes of 
a few hundred kyr (depending on clump masses/densities). Strong 
neutral gas cooling was found to enhance this formation mecha- 
nism, producing denser and longer lived pillar-like structures com- 
pared to models with weak cooling. 

Models such as these for the formation of bright-rimmed 
clouds, gl obules and pi l lars have been considered for many 
vears (e.g. |Pottaschl[l958l : iMarshI Il97d : IBodenheimer et al.ll 19791: 
Sandford et al.lll982h : much of this work is summarised by lYorkd 



J 98^. The RDI of a photoionised clump and its subseq uent 
acceleration and evolution was calculated analytically (.Bertoldil 
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Il989l: iBertoldi & McKed [ 199(|) and subsequently numerically 
bv iLefloch & LazareffI (Il994h . IWilliams et alJ (l200lh considered 
a range of axisymmetric models showing that multiple scenar- 
ios could form long-lived pillar-l i ke structures. It ha s been shown 
dKessel-Devnet & Burkerj l2003l : iMiao et all l2006l : IPound et all 
I2OO7I : iBisbas etal.ll201Qh that RDI of single clumps can gener- 
ate cometary globules and trigger gravitational collapse, but it 
is more difficult to form pillars like those in Ml 6 because gas 
must accumulate to a high density in the shadowed tail region. 
iLim & Mellemal (l2003h showed how photoionisation of multiple 
clumps which partially shadow each other leads to de nse gas accu- 
mulating in shadowed tail regions; it was suggested bv lPound et al.l 
(I2OO7) that multiple clumps are required to quant itatively match 
observations of the pillars in Ml 6. Recent models ( Mellema et al.l 



2006bl : lRaga et al.ll2009l : rLora et al.ll200^ : lGritschneder et al.ll20o'9r 
2010, ML 10) showed that elephant trunks can form quite naturally 



from the photoionis ation of a clumpy densi ty field under a range 
of initial conditions. iGritschneder et aP (I2OIQ) extended our previ- 
ous work (ML 10) which used static initial conditions by consid- 
ering the photoionisation of dynamic density fields generated by 
isothermal decaying turbulence, measuring the formation of pillar- 
like structures as a function of turbulent Mach number. 

Filamentary structure in an apparently helical geometry was 
found in the shadowed tail regions of a number of elephant trunks 
dCarlqvist et al.lll998Ll200l [2003). and it was suggested this could 
arise due to twisting of magnetic field lines. Velocity profiles for 
som e trunks were sho wn to be consistent with solid-body rota- 
tion (iGahm et al again consistent with a magnetic origin of 
the observed structure, although the actual magnetic field orienta- 
tion and strength has not yet been measured for these structures. 
The magnetic field in cometary globules has been measured by 
optica l polarimetry dSridharan et al.ll 19961 : iBhatdl 19991 : iBhatt et all 
I2OO4I) . showing in two cases a field orientation along the head-tail 
axis of the globul e, but in one case Ighatt 1999) a perpendicular 
field was found. ISugitani et al.l (l2007h used near-IR polarimetry of 
background stars to measure the magnetic field in Ml 6, finding an 
ordered large-scale field in the H II region, but within the pillars 
the field is aligned with the pillar axes and significantly misaligned 
with the ambient field by ^ ^ 30-40°. They suggest this should 
constrain the magnetic field strength since it has not been strong 
enough to resist reorientation during the formation and evolution 
of the pillars, a suggestion we explore in more detail in this work. 

Theoretical calculations of the effects of magnetic fields on the 
expansion of H II regions were first considered b y Lasker (1966b). 
Using analytic calculations iRedman et al.l(ll998h studied ionisation 
fronts with a plane-parallel magnetic field in the plane of the ionisa- 
tion front. They found that the velocity separation between R-type 
and D-type solutions decreases with increasing field strength, and 
that D-critical ionisation fronts also advan ce into the neutral gas 
more rapidly for increasing field strength. In I Williams et al.l (l2000h 
jump conditions for ionisation fronts with oblique magnetic fields 
were presented, together with ID numerical models showing how 
an ionisation front could progress from fast-R-type through fast-D- 
type, slow-R-type, and finally to slow-D-type. These extra modes 
are allowed because fast and slow shocks detach from the ionisa- 
tion front at different times and propagate into the neutral medium. 
This work was extended by Williams &Dvson ( 2001), who cal- 
culated the internal structures of stationary ID ionisation fronts. 
They found that oblique fields could produce significant transverse 
velocities and regions where the transverse field component is sig- 
nific antly mod ifi ed. 

IWilHamsl (l2007h used 2D slab-synmietric radiation- 



magnetohydrodynamics (R-MHD) simulations to study the 
photoevaporation flows from magnetised globules. In these models 
clumps with an initial density of nn = 2 x 10^ cm~^ were placed 
in an ambient medium 10 x less dense. For simplicity a uniform 
magnetic field with various strengths and orientations was used. 
Plane-parallel radiation was assumed, with a thermal model where 
the temperature relaxed to a value between 100 K and 10 000 K 
according to its ionisation fraction. It was found that for a weak 
field the ionised gas pressure dominates the dynamics and the field 
was swept into a configuration where it was parallel to the column 
of neutral gas behind the dense clump. For a sufficiently strong 
field, however, the field determined the dynamics and made the 
flow almost one-dimensional along field lines. 

The first 3D R-MHD calculation including non- 
equilibrium^^ transfer was performed 
by iKrumholz. Stone & Gardinej (l2007h . They used the Athena 
code with a new ray-tracing module to simulate the expansion of 
an H II region around a point source in a uniform magnetic field. 
The overall exp ansion is now ax isymmetric rather than spherically 
symmetric (cf. iLasked Il966ah with a dense shell forming in 
directions parallel to the field, and a possibly numerical instability 
developing for expansion perpendic ular to t he field. 

Using similar methods, He nnev et al.l (1200 9) model the pho- 
toionisation of a dense clump of gas in 3D with an initially uniform 
magnetic field. They use the photon-conserving C^-ray ray-tracing 
algorithm (Mellema et al. 2006a) which allows large timesteps to 
be taken without loss of accuracy. They found that the evolution of 
a photoionised globule can be significantly altered by the presence 
of a strong field, and in some cases a recombining shell forms at the 
termination shock of the photoevaporation flow when it is confined 
by a transverse field. The implosion phase is strongly asymmet- 
ric since the clump compresses much more readily along field lines 
than across them. These authors introduce a detailed thermal model 
to model the dynamics as realistically as possible for conditions in 
the Orion Nebula, finding that X-ray heating keeps the neutral gas 
temperature at T > 50 K, but the cooling in neutral gas is some- 
what stronger than that considered in ML 10. 

In this work we add magnetic fields of various strengths and 
orientations to some of the models considered in ML 10 to study 
their effects on the dynamics of the dense neutral gas. The sim- 
ulation code is described in ^ ^2. 1 1 describes the MHD dynam- 
ics algorithm; ^2.2| reviews the microphysical processes included; 
^2. 3 1 presents a brief comparison of results obtained with our code 
to the test problem described by Krumholz et al. (2007). The 3D 
simulations presented here are introduced in |3]and our results pre- 
sented in SI ^4.11 shows the evolution of the projected gas den- 
sity and magnetic field orientation; ^4.21 shows the emission from 
ionised gas; ^4.31 evaluates boundary effects in strong field simu- 
lations. A calculation explaining the presence of a bright ionised 
linear ridge/ribbon in the strong field simulations is presented in 
|5] These results are compared to observations of H II regions and 
their magnetic fields in |6] where we also discuss our work in the 
context of other recent computational research. Our conclusions are 
presented in |7l and some technical details and tests for the simu- 
lations are given in the appendices. 



2 NUMERICAL METHODS 
2.1 MHD algorithms and tests 

These calculations are performed on a uniform finite-volume grid, 
with the hydrodynamics and ray-tracing/microphysics modules as 



© 0000 RAS, MNRAS 000, 000-000 



Magnetised Pillars and Globules 3 



previously described in ML 10. The simulations presented here 
use the same code, but with an MHD module coupled to the 
ray- tra cing. The integration scheme was described by Fall e et all 
(Il998h and is second order accurate in time and space, using the 
symmet ric van Albada slope li miter to ensure monotonicity near 
shocks dvan Albada et To calculate the MHD fluxes a 

number of Riemann solvers were trie d; we use the Roe solver in 
conserved variables described by Car go & Gallicd (Il997l) as it was 
the most robust for these ca lculations (with co rrections to typos ob- 
tained by comparison with lStone et al.|[2008b . Magnetic field units 
are used in the code such t hat factors of Atv d o not appear in the dy- 
namical equations (see e.g. Falle et al. 1998); c.g.s. units can be ob- 
tained by scaling the field strength by a/Itt. To avoid confusion we 
will always quote field strength in c.g.s. units (usually micro-Gauss, 
/iG) rather than the code units. We use notation such that the total 
energy per unit volume is E — ^pv"^ + Pg/{l — 1) + Pm, being 
the sum of the kinetic, internal (thermal), and magnetic energies. 
The gas pressure is pg, magnetic pressure is pm = /Stt, and the 
ideal gas adiabatic index 7 is set to 5/3. Multid imensional calcu- 
lations use the artificial diffusion prescription of Falle et all (Il998l) 
to fix the 'carbuncle problem' and to prevent rarefaction shocks de- 
veloping; typically a viscosity parameter 77^ = 0.1 is used. 

The code uses cell-centred values for vector quantities; the 
divergence constraint in the magnetic field i s addressed using th e 
'mixed-GLM' divergence cleaning method of Dedner et 
This algorithm introduces an extra scalar field which couples to 
V • B, advecting and diffusing divergence errors so that they do not 
build up near shocks or other discontinuities. While it does allow 
non-zero divergence of B (which can introduce small force errors), 
the method has some advantages. (1) It is fully conservative in all of 
the physical variables; only the evolution equation for the unphysi- 
cal scalar field, ip, contains a source term. (2) The total energy and 
magnetic field updates are calculated in the same algorithmic step, 
ensuring they are consistent. An energy correction is still required, 
but this can be calculated from the fluxes generated by divergence 
cleaning. This means the internal energy (as a fraction of the to- 
tal energy) is treated quite accurately, making the code robust for a 
wide range of field strengths. (3) It involves relatively little compu- 
tational and memory overhead. 

A number of ID and multidimensional adiabatic MHD 
code tests have been performed to validate the code and 
for comparison to other astrophysical MHD codes. Re- 
sults from a number of these test problems are shown 
in Appendix |a1 and can also be found in more detail at 
|http : //www .astro . uni-bonn . de/~ jmackey/ jmac/| 
together with HD and photoionisation test problems (some of 
which were presented in ML 10). 

Three modifications have been made to the mixed-GLM diver- 
gence cleaning method. (1) The parameter Cr is set to Cr = 4 Ax 
(where the grid cell dia meter is Ax) inste ad of the recommended 
value of Cr = 0.18 (see lDedner et al IE002I . eq. 45). A typo in their 
paper (A. Dedner, private communication) meant Cr should have 
been defined as Cr = Ax/0.18, similar to the value used here. 
(2) The algorithm induces extra magnetic field transport across cell 
boundaries; the energy associated with this is accounted for by an 
extra term in the energy flux: 

F[E]^ F[E] + F[B^]B^,m . (1) 

Here F[X] represents the flux of conserved variable X across a 
cell boundary, Bx is the component of B normal to the cell bound- 
ary, and Bx,m is the value of Bx in the interface state calcu- 
lated using the mixed-GLM algorithm CDedner et al...200Z . eq. 42). 



This correction removed gas pressure dips which appeared ahead 
of oblique shocks in axisymmetric models of magnetised jets. (3) 
Motivated by numerical difficulties encountered for subsonic flow 
across boundaries in magnetically dominated media, the boundary 
condition for the scalar field ip in the mixed-GLM method is modi- 
fied for zero-gradient boundaries. This is described in Appendix iBl 
showing the improvement obtained with the new boundary condi- 
tion. 



2.2 Coupling MHD to photoionisation 

We briefly review the microphysics algorithm here; the reader is 
referred to ML 10 for further details. The dynamics and micro- 
physics are solved separately in turn by operator splitting: in a 
given timestep the adiabatic MHD equations are solved and up- 
dated first (with the ion fraction advected with the gas flow), 
and then the internal energy per unit mass (cint) and ion frac- 
tion of Hydrogen (x) are updated by an adaptive sub- stepping 
integration to a user- specified relative accuracy. In the absence 
of photoionisation this uses a 5th order Runge-Kutta integra- 
tion. Ray-tracing uses the short characteristics method with the 
photon-conserving discretisation of the photoionisation rate given 
bv lMellema et al.l (l2006aD. When a cell is strongly photoionised ex- 
plicit integration is very unstable and so the 'constant elect ron den- 
sity' approxima tion of | Schmidt- Voigt & KoeppenI (Il987h is used 
(see also Mellema et al.l l2006al) to analvticallv integrate the ion 
fraction when x > 0.95. The time- averaged fraction of photons 
which pass through the cell is calculated simultaneously to obtain 
the appropriate cell optical depth for calculating rays to more dis- 
tant grid cells. The algorithm used in ML 10 was designed to work 
in an identical way for R-MHD as for R-HD. As such, the non- 
dynamical test problems presented in ML 10 produce identical re- 
sults when run with a non-zero magnetic field. 

We continue to use the on-the-spot approximation for diffuse 
radiation; this remains a potential limitation of our results which 
should be checked with future work. Recent calculations including 
diffuse radiation ( Raga et al. 2009; Williams & Henney 2009) have 
shown that its role in the dynamics of H II reg i ons m ay not be as 
significant as was suggested by e.g. Ritzervel dl (l2005h . We plan to 
test this in d ynamical sim ulations using an algorithm such as that 
described bv iKuiper e t al. (2010); it should be possible to extend 
their method to treat ionising photons. 

Radiative cooling processes are modelled using the C2 cool- 
ing function from ML 10 - in photoionised gas the coolants con- 
sidered are collisionally excited emission from Oxyge n and Ni- 
troge n dOsterbrocS Il989h and recombining Hydrogen (iHummej 
I1994I) . and neutral gas cooling is assumed to follow Newton's Law 
of exponential cooling with a cooling time tc = lOkyr. Photoion- 
isation assumes a monochromatic ionising source with 5.0 eV of 
thermal energy added to the gas per ionisation. Collisional ionisa- 
tion is calculated using the rate of lVoronovl ( Il997h . 

2.3 H II region expansion with a magnetic field 

iKrumholz et al 1 (I2OOI hereafter KSG07) added a photoionisation 
module to the Athena code to calculate the 3D R-MHD expansion 
of an H II region into a uniform medium with an initially uniform 
magnetic field. We have modelled this problem both in axisymme- 
try and in 3D to compare results from our code to those of KSG07. 
Our 3D model used a smaller domain of diameter 14 pc resolved by 
160^ cells (12 per cent larger cell diameters than KSG07). The only 
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Config. 


Object 


X 


y 


z 


So 


ro 


M 


1,2 


^min 


1.5 


-1.5 


-1.5 


n/a 


n/a 


n/a 


1,2 


^max 


6.0 


1.5 


1.5 


n/a 


n/a 


n/a 


1,2 


Source 











n/a 


n/a 


n/a 


1 


CI 1 


2.30 








500 


0.09 


28.4 


1 


CI 2 


2.75 





0.12 


500 


0.09 


28.4 
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CI 3 


3.20 





-0.12 


500 


0.09 


28.4 


2 


CI 1 


2.30 





-0.174 


250 


0.12 


33.7 


2 


CI 2 


2.30 





0.174 


250 


0.12 


33.7 


2 


CI 3 


2.73 








500 


0.09 


28.4 


2 


CI 4 


4.20 


0.12 


0.12 


400 


0.12 


53.9 



Table 1. Source and clump properties in R-MHD simulations. Configura- 
tion 1 is model 17 from ML 10 and configuration 2 is a modified version 
of model 18 with an extra clump added further from the source. The do- 
main minimum/maximum coordinates are measured in parsecs relative to 
the source, which is at the centre of the y-z domain. Clump overdensities, 
6o, are relative to the background number density of nu = 200 cm~^; 
their Gaussian scale radii, ro, are in parsecs and their total masses, M, in 
units of M0. 



Figure 1. R-MHD model of 3D H II region expansion. Results are shown 
at 1.58 Myr for a slice through the symmetry axis of the 3D simulation. 
Coordinate axes show position in parsecs, and gas number density is shown 
on a logarithmic scale from nu = 1-10^ cm~^ as indicated. Velocity 
vectors are also shown where vector size is linearly proportional to velocity, 
V. The maximum velocity is ^;max = 6.7 km in this figure. A white 
contour marking ion fraction x = 0.5 is also shown. The 14.2 fiG magnetic 
field is initially horizontal (parallel to the symmetry axis). 



other difference in our simulation was to add a low level of noise 
to the otherwise uniform initial density and gas pressure fields. An 
initial ambient temperature of T = 100 K was used rather than 
the 11 K used by KSG07 and we used the C2 cooling function 
described above, but tests showed this does not affect the results 
significantly. 

Gas density and velocity field are shown for a slice through 
the x-y plane at t ^ 1.58 Myr in Fig. [T] Comparison to fig. 18 
of KSG07 shows that very similar results are found here, notably 
the dense shell which forms for expansion in the =bx directions and 
its absence in perpendicular directions. We also find the instabili- 
ties which were noted by KSG07 in the perpendicular direction; we 
suggest they are real (i.e. not numerical artefacts) and related to the 
R-HD ionisation front instabilities studied by previous authors (e.g. 
iGarcia-Segura & Francol[l996l : IWilliamsll2002l : IWhalen & NormanI 
12008*) on the basis of tests which showed the instabilities are only 
present when neutral gas can cool strongly. The velocity arrows 
show that the magnetic field significantly alters the ionised gas dy- 
namics: at late times the strong outflow along field lines transforms 
the ionisation front to a recombination front in the dzx directions. 
We plan to study these features in more detail in future work. 



3 3D R-MHD SIMULATIONS 

To focus specifically on the effects of magnetic fields on the pho- 
toionisation process, the dense clump configurations in models 17 
and 18 from ML 10 were used rather than the random clump dis- 
tributions in models 1-16. The two clump configurations as well 
as the simulation domain and radiation source properties are de- 



scribed in Table [T] Configuration 1 is basically the same as model 
17 in ML 10; configuration 2 is changed slightly - the two front 
clumps are slightly larger and further apart, and an extra clump 
is added further from the source. The radiation source located 
at the origin has the same properties as before: an ionising pho- 
ton luminosity = 2 x 10^° with monochromatic photon 
energy hvo — 13.6 eV = 5.0 eV. The simulation domain used 
here is twice as large in every dimension as was used for mod- 
els 17 and 18 in MLIO; this is because boundary effects were 
found to be niuch m ore significant with strong magnetic fields (cf. 
iHennev et al.ll2009h . The simulation domain contains 384 x 256^ 
cells, giving a physical resolution of Ax 0.012 pc per cell. Re- 
sults of 3D R-MHD simulations with various domain sizes are de- 
scribed in Appendix [eland demonstrate that the domain used here 
is sufficient to prevent significant boundary effects for all but the 
strongest fields modelled. 

These models were run with zero-gradient boundary condi- 
tions, which perfectly model supersonic outflow and damp reflected 
waves in subsonic outflow. There is nothing to stop inflow, however, 
if that is what the solution tends towards, and some of the R-MHD 
simulations do develop strong inflows at later times. An alterna- 
tive boundary conditi on was imposed in the 2D slab- symmetric R- 
MHD simulations of Williamsl J2007h : when the velocity was out- 
ward the usual zero-gradient boundary was applied, but if the veloc- 
ity changed sign then inflow was suppressed. This 'only-outflow' 
boundary condition was also tested in our simulations to study how 
significantly inflows affected the solution, and in particular its ob- 
servable properties. This was initially tested using axisymmetric 
simulations with a magnetic field parallel to the radiation propa- 
gation vector; the results showed that for a weak or medium field 
up to B ^ 50 /iG the inflow made little difference. Its only ef- 
fect was to confine the photoevaporation flow from the clump to a 
smaller volume by the inflow's ram pressure. The higher mean den- 
sity of ionised gas (and associated stronger recombination) had no 
discernible dynamical effect. With a strong field of B ~ 160 /jG, 
however, the inflow could overrun the photoevaporation flow and 
the ionisation front, transforming it into a recombination front and 
hence strongly affect the solution. For this reason the strong field 
simulation R8 described below was run with both a zero-gradient 
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(as R8) and an only-outflow (as R8a) boundary condition in the 
direction back towards the radiation source. 

3.1 Description of simulations 

The simulations are listed in Table [21 Rl and Rll are R-HD sim- 
ulations of clump configurations 1 and 2 respectively, while the 
other simulations are R-MHD calculations of the same clump con- 
figurations with varying field strengths and orientations. The field 
strengths were approximately B 18, 53, and 160 /xG and are 
referred to as weak, medium, and strong respectively. R2-R10 are 
R-MHD simulations of configuration 1 : R2-R4 are weak field mod- 
els, R5-R7 have a medium field, and R8-R10 a strong field. R12- 
R15 are R-MHD simulations of configuration 2: R11-R14 have a 
weak field and R15 a medium field. 

For the field orientation, 'parallel' denotes parallel to the ra- 
diation propagation vector along the centre of the simulation (B = 
5ox), and 'perpendicular' denotes a field in the y-z plane. Field 
orientations were chosen to be either parallel or perpendicular, or 
oriented 80° from the x-axis and 50° from the y-SLxis in the y-z 
plane. The simulations with the latter field configuration (R3, R4, 
R7, RIO) produced very similar results to the perpendicular models 
(R2, R5, and R8) for configuration 1 so they were not run for as 
long as the other models and were not run at all for configuration 2. 

Whether a field is weak or strong is largely determined by its 
dynamical importance, set by the plasma parameter (3 = STrpg/B"^, 
the ratio of thermal to magnetic pressure. Since ionised gas in 
H II regions is largely isothermal, as is dense molecular gas, the 
thermal pressure is approximately proportional to density. Hence 
a field of 50 /xG could be strong or weak depending on whether 
the gas is ionised or neutral, and on its density. The plasma pa- 
rameter is shown in Table [3] for a range of gas pressures encoun- 
tered in the photoionisation simulations: (1) the initial conditions 
have a constant pressure of pg — 1.38 x 10~^^ dynecm"^ (or 
Pg/k^ — 10^cm~^K); (2) gas at the background density of 
riH = 200 cm~^ and the ionised gas temperature of T ~ 8 000 K 
has Pg — 4.42 x 10~^° dyne cm~^, and (3) the peak pressure typ- 
ically encountered in the simulations is pg ^ 10~^ dyne cm~^. 
For weak field simulations the gas pressure clearly dominates and 
for the strong field the situation is reversed; for the medium field 
case the initial conditions are magnetically dominated but, once 
ionised, the gas pressure is larger. Thus it is expected that the weak 
field results will largely follow the R-HD results, the strong field 
simulations should show very different behaviour, and the medium 
field models will lie somewhere in between. The gas pressures in 
these simulations are comparable to those estimated for the pillars 
and their environment in M16 (see the discussion in MLIO), but 
it should be borne in mind that other massive star-forming regions 
can have significantly higher (or lower) gas pressures. 



4 RESULTS 

Three main observable consequences of th e magnetic field are ex- 
pected (IWilliamsl2007l : lHennev et al.l2009L cf.): (1) For weak mag- 
netic fields the field orientation will be be changed by the dynam- 
ics of the photoionisation process. (2) When the field is sufficiently 
strong the density structure of the neutral gas will be significantly 
altered because gas can only move along field lines; RDI produces 
a sheet rather than an axisymmetric compression. (3) Strong mag- 
netic fields will confine the photoevaporation flow changing its ob- 
servable properties. We first plot projections through the simulation 
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Table 2. Clump and initial magnetic field configurations, boundary condi- 
tions applied, and final simulation times for the R-MHD simulations. The 
two clump configurations indicated in column 2 are described in Table [T] 
For the boundary conditions (BCs, column 3), 1 refers to zero-gradient, 
2 refers to only-outflow, and 3 is a hybrid where all boundaries are zero- 
gradient except for the boundary facing the source which is only-outflow. 
Field strengths are given in units of /xGauss. Some models have slightly 
off-axis field orientations to prevent any potential grid-alignment numerical 
effects. The simulations were evolved to the time shown by tgim (in kyr). 
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Table 3. Plasma beta {(5 = 87rp^/B^) values for gas pressures encountered 
in R-MHD simulations, with /3 <C 1 indicating magnetically dominated re- 
gions. The field strengths quoted are the three initial field strengths applied 
in the simulations. Gas pressures (dyne cm- ^) are representative of the 
initial conditions (2nd column), ionised gas at the background density (3rd 
column), and the typical maximum pressure near the ionisation front (4th 
column). 

domain showing column density of neutral gas and the projected 
magnetic field. Following this we show emission maps in recom- 
bination radiation. R2, R5, and R8 proved most useful for showing 
the effects of increasing field strength so we focus most of our anal- 
ysis on these three simulations. 

4.1 Projected density and magnetic field 

The column density is calculated as in MLIO by integrating the 
neutral gas number density along the line of sight (LOS). Note that 
because we do not consider molecules explicitly, the column den- 
sities should be divided by 2 to give N(H2). The projected mag- 
netic field is more difficult to calculate since it must be a weighted 
integral, for example to calculate the polarisation of background 
starlig ht indu ced by aligned dust grains (cf. observations of Ml 6 
bv lSugitani et allboovh . For our integration we assume a constant 
gas-to-dust ratio and weight the integral by gas density. To allow 
for the possibilit y that grain alignme nt may be less effective at high 
densities (e.g. Goodm an et al .Il995h we limit the density weighting 
to a maximum density, rimax = 2.5 x 10^ cm~^; this is rather ad- 
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Figure 2. Projections along the y-axis for simulations R2 (top), R5 (centre), and R8 (bottom) at simulation times 100 kyr (left) and 200 kyr (right). Only part 
of the simulation domain is shown; positions are shown in parsecs relative to the source. Neutral gas column density is plotted on a logarithmic scale (cm~^) 
as indicated. Projected magnetic field orientation (Eq.[2} is indicated by the black lines. 



hoc but it has a very limited effect on the projected field orientation 
as long as rimax ^ 10^ cm~^. The projected field may be calcu- 
lated by integrating the Stokes Q and U parameters along the LOS 
and subsequently r ecovering th e field orientation using trigonomet- 
ric relations (see Arthur et al. 2010). With the approximations de- 
scribed above the projected Q and U parameters of the magnetic 
field are given by 



min[nH(s),nn 
(U) = / min[nH(s),nn 

J s^O 



Bl 



2B1B2 



ds , 



(2) 



where s is the LOS and the field perpendicular to the LOS is 
Bp = [Bl, B2]. Technically this integration has units of G.cm~^ 
but the normalisation was not considered because we only use the 
projected orientation of the field in this work and not its magni- 
tude. In principle the relationship between polarisation and mag- 
netic field is a very complex function of temperature, density, local 



radiation field and dust properties fe.g. lSpitzeJ 19981 ) but the micro- 
physical processes included in our simulations are not sufficient to 
model this in any detail. 

Results from the perpendicular field models R2, R5 and R8 
are shown in Figs. El and [3] for simulation times 100, 200, 300 
and 400 kyr. The LOS is the simulation y-axis so the initial field 
in all three simulations is almost vertical and perpendicular to the 
LOS. At 100 kyr the leading dense clump is just past the point of 
maximum compression due to RDI. The partially shadowed clumps 
are being compressed asymmetrically and are not yet at maximum 
compression. At this stage the structure of the dense gas is very 
similar in all three models, although it should be noted that the 
low density tail of the strong field model (R8) is actually a thin 
sheet seen edge-on, whereas it is closer to axisymmetric in R2. Al- 
ready the magnetic field has been significantly altered in R2 due to 
RDI; this is less apparent in R5 while the field in R8 is almost un- 
affected. By 200 kyr the clumps have re-expanded somewhat (lim- 
ited because much of the compressional heating was radiated away) 
but the rocket affect has taken hold and the first clump has almost 



© 0000 RAS, MNRAS 000, 000-000 



Magnetised Pillars and Globules 7 



3.0e+19 



3.0e+20 



3.0e+21 



3.0e+22 



3.0e+23 




0.0 



3.0 4.0 3.0 4^0 5.0 

3.0 4.0 




0. 



mmm>% 



'■mmiMim 

3.0 4.0 

Figure 3. As Fig. [2 




but for simulation times 300 kyr and 400 kyr. 



fully merged with the second. At this stage the magnetic field re- 
orientation in R2 is more pronounced and in R5 is beginning to 
become significant. There are also differences in the morphology 
between the three models, quite slight between R2 and R5 but sig- 
nificant for R8. After 300 kyr the differences have become even 
more pronounced although the three clumps have basically merged 
to one in all three simulations. By 400 kyr R2 and R5 have become 
cometary globules rather than pillar-like structures. R5 is also frag- 
menting - the small protrusion at the top of the dense region even- 
tually detaches from the main clump and is rapidly accelerated. R8 
at this stage looks completely different to the other models. Low 
density gas has recombined behind the very broad ionisation front 
formed by gas expansion along the field lines. The field remains 
almost unchanged from its initial value even at this late stage. 

The field evolution is shown in a more quantitative way in 
Fig.lH where the ratio of the mean (volume averaged) parallel field 
(|Bcc|) to perpendicular field (y^S|~T~Bf) is plotted as a func- 
tion of time for both the full simulation domain and for only those 
cells with gas density nn ^ 5 000cm~^. This ratio changes rel- 
atively little in the full domain because ionised gas takes up the 
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overwhelming majority of the simulation volume. For the dense 
gas, however, a very strong effect is seen in the evolution as the 
field strength increases. The weak field in model R2 is swept into 
alignment with the pillar by the dynamics of RDI and the rocket 
effect, but a sufficiently strong field (here ~ 160 /xG in R8) is un- 
affected. 

In the parallel field models R6 and R9 there was very little 
field evolution because the rocket effect acts in the same direction 
as the initial field orien tation. The re s ults ob tained are comparable 
to the SOOL model in iHennev et al.l (l2009h . except that for R9 a 
strong inflow from the boundary closest to the source overran the 
ionisation front after t ^ 200 kyr, leading to a stable ID recombi- 
nation front. This would n ot be expected to form if an only-outflow 
boundary is imposed, as in IWilliamsl ( [2007h . although we note that 
for ID photoionisation where ionised gas is confined to stay in the 
LOS back towards the source, a recombination front is a more sta- 
ble and natural solution than an ionisation front. 

4.2 Emission maps in recombination lines 

The gas emission in recombination radiation (e.g. Ha) was also 
calculated from the sin gulation outputs as in ML 10, using the 
e missivity formu l a froni iHennev et alj (l2005h and dust opacity as 
in lMellema et al.l (l2006bl) to attenuate radiation along a LOS. Emis- 
sion maps are shown at times 100, 250, and 400 kyr for models 
R2, R5, and R8 in Figs. [5l[6l and [T] Here, as in the column density 
maps, there are small differences between R2 and R5, but R8 shows 
very different emission. This can be understood from Table[3]which 
shows that only R8 is magnetically dominated in ionised gas. The 
initial field is along the LOS for the left-hand projections (LOS is 
simulation z-axis) and vertical for the right-hand projections (LOS 
is simulation ^-axis). 

The field is so weak in R2 that the evolution is very similar 
to the purely hydrodynamical evolution in Rl (cf. ML 10) except 
that the low density shadowed tail region is not quite axisymmet- 
ric. A strong photoevaporation flow from the front clump expands 
spherically with velocity v 30-35 km until the ram pressure 
equals the total pressure of the ambient medium at which point a 
standing shock is established. For the full evolution of this model 
the dense gas at the ionisation front is by far the brightest structure 
in recombination emission. There are small differences in model 
R5, notably that the photoevaporation flow is confined to a smaller 
volume by the higher ambient total pressure, but the evolution is 
largely the same. 

In R8 the total pressure of the ambient medium is magnetically 
dominated and much higher than in R2 or R5. This confines the 
photoevaporation flow much more effectively so that the region of 
spherical expansion is significantly smaller. Instead shocked pho- 
toionised gas is diverted along field lines into a bar- shaped lin- 
ear structure with a density 2.5-4 x the background density. In 
IHennev et al this structure is termed a dense ribbon to avoid 

confusio n with photon -dominated regions such as the Orion Bar 
(see e.g. lO'Dellll200lh . We will avoid the term 'bar' for the same 
reason, instead refe rring to the structure as a ribbon or ridge (cf. 
iBohigas et al .1120041) . This dense ionised gas is quite spatially ex- 
tended and so has a comparable recombination radiation intensity 
to the bright rim at the ionisation front. Because of the high col- 
umn density in the ribbon there is significant recombination in its 
shadow at later times. This evolution is si milar to that seen in the 
strong field models of IHennev et al.l (l2009h . 

The left-hand plots for R8 (with B along the LOS) show 
wave-like features which appear to be partially developed Kelvin- 



Helmholz instabilities. These move back towards the midplane 
(y = 0) due to magnetic tension and have collided by t = 400 kyr 
(Fig. |7]); similai^but more pronounced evolution was found by 
IHennev et al.l (l2009h in their simulations. 

In simulation R8 at 400 kyr the dense ridge/ribbon has become 
optically thick, enhancing recombination in its shadow and dramat- 
ically weakening the photoevaporation flow. The result of this is 
that the standoff distance between the ridge and the original ioni- 
sation front of the dense clump has decreased almos t to zero. This 
was found to be a short-lived phase in the models of IHennev et al.l 
(2009) and we expect the same to be true here because a strong 
photoevaporation flow is required to generate the overdense ionised 
gas. 

4.3 Boundary effects for simulation R8 

It was found that a reasonably strong inflow of gas developed at the 
boundary x = 1.5pc in R8 at late times, ranging from Vx = 1- 
3.5 km s~^. The ram pressure from this inflow could raise the den- 
sity in the shocked region and hence affect the resulting evolu- 
tion, so the simulation was repeated with this boundary set to only- 
outflow, thereby preventing the inflow developing. This simulation, 
denoted R8a, is shown in Fig. [8] at times 100, 250 and 400 kyr, 
again in recombination line emission. Comparison to Figs. |5] and [6l 
shows that the evolution is almost identical up to t = 250 kyr ex- 
cept that the bright ridge is slightly further from the ionisation front 
because now there is no ram pressure confinement. Ait — 400 kyr 
the differences are more significant (see Fig. [3: the ionised ridge 
is broader and less well-defined in R8a than in R8, and while the 
dense neutral gas has a similar structure there are small differences. 
The general features remain the same, however, and the confined 
photoevaporation flow still shields the ionisation front significantly 
and outshines it in recombination radiation. 

A similar but weaker inflow was also found in simulation R5. 
The only effect this has on the results is that the ionised ambient 
medium is denser than in R2, as seen by comparing the mean emis- 
sion between the two simulations in Figs.[5HZl 



5 IONISED RIDGE FORMATION 

The presence of a bright ionised ridge/ribbon of dense gas in sim- 
ulation R8 which is absent in Rl, R2, and R5 is easily understood 
quantitatively by evaluating the jump conditions for the isothermal 
termination/reverse shock of the photoevaporation flow (the tem- 
perature is always in the range 7500-8500 K in the ionised gas). 
If we assume the system has reached an equilibrium, then the to- 
tal pressure in the upstream (I), shocked (II), and downstream (III) 
regions will be equal. To avoid dealing with the spherical expan- 
sion of the photoevaporation flow we consider region I to consist 
only of the conditions immediately upstream from the shock (re- 
ferred to below with subscript '1'). If we further assume that the 
only significant magnetic field is transverse to the velocity of ex- 
pansion we can consider only its pressure effects and define a speed 
= / (Sttp) — pm/ p so that the total pressure in a region is 

Ptot = Pram + Pg + Pm ^ p{v^ + + 6^) (3) 

being the sum of the ram, thermal and magnetic pressures respec- 
tively. 

We wish to calculate the density in region II as a function 
only of the upstream velocity vi and the ambient medium prop- 
erties (63, ^^3). Only when the ratio P2/ p?> > 1 will a bright ribbon 
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Figure 5. Projected recombination radiation maps of simulations R2 (top), R5 (centre) and R8 (bottom) for projections along the z-axis (left) and y-axis (right), 
here att = 100 kyr. The initial field is along the z-axis for all three simulations. Energy flux on the indicated logarithmic scale has an arbitrary normalisation. 
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Figure 6. As Fig.[5]but at t = 250 kyr. 
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or ridge be observable. Thu s the equations to be solved are (cf. 
Ide Hoffmann & TelleJ 19501) 

kiiv^, + + bl) = k2ivl + + bl) = (vl + + bl) , (4) 

where ki = pi/ps and /c2 = p2/p3- The isothermal sound speed 
is = Pg/ p — kT I (pnip), where p is the mean mass per particle 
in units of the proton mass rrip and is p = 0.5 in ionised gas in 
our models. For T ^ 8300 K this gives c = 11.7kms~\ In all 
simulations the magnetic field is negligible in region I (bi <^ c) and 
smaller than the other terms in region II (62 < c). The flow velocity 
in the reverse shock reference frame is vi 33-35 km s"^ in all 
simulations giving a Mach number M = vi/c 3. 

The hydrodynamic jump conditions for the isothermal reverse 



shock are 



p2 



V2 = c/M 



(5) 



When 61 but 62 > the density jump is given by the solution 
to the quadratic equation 



P2 



2 I 2 



+ 



vl 



pij c^ + bl c^ + bl 



: 



(6) 



This equation has real solutions only for some values of 62, leading 
to the constraint 



M2 / 



(7) 
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Figure 7. As Fig.[5]but at t = 400 kyr. 
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Figure 8. Similar to Fig.|5] but showing simulation R8a at t = 100 kyr (top), t = 250 kyr (middle) and t = 400 kyr (bottom). The left-hand plots show 
projection along the z-axis and the right-hand plots along the y-axis. These can be compared to R2, R5, and R8 in Figs.|5]|6]and|7] 



For M = 3 this gives 6| ^ 16c^/9 but in the simulations it is 
actually significantly less than this limiting value. 

Using these equations it is easy to show that in the hydrody- 
namic limit /cs = M^/{1 + M^) 0.9 and ki = 1/(1 + M^) - 
10. No bright overdensity is formed, and the photoevaporation flow 
shocks at a density about 10 x below the ambient density. This is 
indeed seen in the R-HD simulation Rl and also the weak perpen- 
dicular field simulation R2. 

With bi — = and 63, 62 > we obtain 



k2 



hl + c' 



1 + {bs/cf 



(c/Mr+c^ + bl- (^)Vl+Mi(l-i.)^ 



k2 ^ 



1 + M2 



1 + 



(8) 



where the inequalities come from the limiting values for 62 given 
by Eq.|7] With M = 3 this gives 



_9_ 

26 



1 + 



^k2^ 



_9_ 
10 



1 + 



(9) 



where = Pg/pm — /bi. The overdensity in the post- shock re- 
gion is therefore determined by the value of ( 1 + 1//3) in the ionised 
ambient medium, showing that the ambient medium must be mag- 
netically dominated to give an over-dense ridge. For R5 /^s = 4.0 
and for R8 /^s = 0.43, leading to maximum overdensities /c2 of 
1.1 and 3.0 respectively. Actual post-shock densities in R5 corre- 
spond very well to this value, and in R8/R8a the ribbon density is 
nn — 550-700 cm~^ depending on the time. The ambient density 
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is riH = 200 cm~^ so the overdensity of k2 = 2.5-3.5 is again 
close to the predicted value. 

Similarly it is easy to show that ki = (1 + l//33)/(l + M^) 
giving ki = 0.125 and 0.33 for R5 and R8. According to this 
equation the photoevaporation flow should shock for simulations 
Rl, R2, R5 and R8 at nn = 20, 21, 25, 67 cm"^ respectively. 
In the simulations at t = 150 kyr the actual densities are nu — 
25, 22, 31, 75 cm"^, whereas at t = 200 kyr they are rm ^ 
18, 19, 30, 100 cm~^, in both cases comparable to the values 
predicted. 

This simple model has obvious limitations, in particular we 
have not specified the form of the boundary between regions II and 
III - it is not a shock, but if the magnetic field strength changes 
significantly then a contact discontinuity forms to equalise the total 
pressure. Also we have not considered any parallel field compo- 
nent which, if present, would alter the pressure balance. Despite 
this the model quantitatively captures the properties of the over- 
dense ridge. Addition of an opposing ram pressure in the ambient 
medium can clearly have a similar effect to a downstream magnetic 
field in that it increases the density at which the photoevaporation 
flow termination shock forms (and hence also the post- shock den- 
sity P2). The differences between R8 and R8a (Figs. [SHE] show that 
even an inflow of V3 — 3kms~^ can have an observable effect. 
Further work will be required to investigate observational differ- 
ences (possibly in the shape and velocity of gas within the ridge 
and/or whether it is sheet-like or bar-like) which could be used to 
distinguish between a magnetically confined bright ridge and one 
which is ram-pressure confined. 



6 DISCUSSION 



This work confirms the suggestion by 'Sugitani et al.l feOO T*) that 
their observations can in principle constrain the field strength in 
Ml 6 when combined with detaile d simulations. On th e evidence 
of these simulations (and those of iHennev et al.|[200^ it appears 
that the magnetic field in the region around the pillars in Ml 6 
cannot be as large as the 160 /xG field in model R8/R8a, and is 
likely |B| < 50 /xG. This can be deduced by comparing observa- 
tions of the field orientation in the H II region and Ha emission 
around t he pillars with the emission maps presented here and in 
iHennev et al. ( 200 ^). The re combination line intensity in the im- 
ages from Hester et al.l (Il996l) clearly decreases with distance from 
the ionisation front and was shown to be consistent with a spher- 
ically expanding photoevaporation flow, which matches only our 
weak and medium field strength models. Additionally no promi- 
nent ribbon/ridge or sheet-like features are present near the pillars. 

Interpreting the the polarisation measurements within the pil- 
lars in Ml 6 as further evidence for a weak field implicitly assumes 
the pillars formed via a mechanism similar to the one considered 
here. Because the rocket effect always accelerates gas away from 
the radiation source and into the shadowed region, we believe that 
this alignment of weak fields with the pillar axis is a generic feature 
of all such m odels, but the initial cond itions for pillar formation are 
not static (cf . iGritschneder et"aDl20 1 Ol) and it is uncertain how well 
correlated the initial magnetic field orientation would be between 
dense and rarefied gas. 

Compared to the simulations of IHennev et al.l (l2009l) . the RDI 
of clumps in our models is weaker leading to less flattened struc- 
tures. The main reason is that here the clumps are considerably 
more centrally concentrated with higher initial densities, making 
them less susceptible to compression. Additionally the clump in 



their model is much closer to the ionising source and subject to 
higher incident radiation flux (only partially offset by the higher 
source luminosity in our model), increasing the effect of RDI. Ther- 
mal physics models and differing geometric effects due to source 
proximity may also have an effect. 

The weak field model W80L in IHennev et"aD (l2009l) should 
be roughly comparable to our model R5. Fig. |4l shows that the ra- 
tio {\Bx\) / {y^B^ + B|) in dense neutral gas remains at ^ 0.3 
for 300 kyr in R5, similar to the ratio in neutral gas shown in 
IHennev et al.l bOOl fig. 12) for W80L. The subsequent increase 
in this ratio in R5 may represent a later evolutionary phase but, be- 
cause of the simulation differences already mentioned, it is difficult 
to quantify this. 

Our results are otherwise in good agreement with those of 
IHennev et al.l (l2009h considering the differences in the simulation 
initial conditions. We also find that that the photoionisation of 
dense clumps is very different in strongly magnetised media than 
in the non-magnetised case, with the dynamics becoming more pla- 
nar than axisymmetric. Photoevaporated ionised gas tends to form 
ribbon-like structures which can (transiently) have significant op- 
tical depths and lead to recombination behind them. Clump com- 
pression along field lines creates flattened sheet-like structures and 
an ionisation front which is far from axisymmetric. 

Despite the simplistic initial conditions used here and in 
ML 10, qualified support for these models comes from the 3D R-HD 
simulations of Gritschneder et al. (2010). They used an isothermal 
self-gravitating decaying turbulence model as the initial conditions, 
varying the turbulent Mach number at which the radiation source is 
switched on. In agreement with the trends we reported they also 
found that to form pillars the density field was required to include 
both large low-density regions and reasonably massive dense re- 
gions of sufficient overdensity to trap the ionisation front. Because 
of the dynamic initial conditions in their simulations they were also 
able to show the dependence on initial gas motions, finding that 
high Mach number turbulence did not produce pillar-like structures 
as successfully as models with Mach number M ^ 4-10 because 
structures formed and dispersed too rapidly. When the Mach num- 
ber was too low a sufficient de nsity contrast to form pillars was 
not obtained. lOritschneder et a l. ( 2010) also found that ionised gas 
pressure is the dominant driver of the dynamics and gravity appears 
to play a smaller role, a n explicit assu mption in our work based on 
previous calculatio ns bvlEsQuivel & R aga (2007). 

In our opinion Gritsc hneder et al.l ((2pi0) overstate the differ- 
ences between our results and theirs - in both models a combina- 
tion of RDI and the rocket effect reinforces pre-existing inhomo- 
geneities in the ISM and generates elongated structure along the 
radiation propagation direction. We do not see our r esults as be- 
ing in conflict with those of Gritschneder et al.' (lOlO); rather that 
they have confirmed many of our conclusions and also extended our 
work by studying pillar formation in a more realistic initial density 
field. These results suggest that the initially static models consid- 
ered in ML 10 and in this work can capture the essential physics 
of the pillar formation and evolution, once the initial conditions of 
dense regions surrounded by a less dense ambient medium have 
been set up. As noted in ML 10, our mo dels can also explain the 
LOS velocities seen in the pillars in Ml 6 ('Pound" 1998V W hite et alj 
1999), although we do not see the helical morphology and appar- 
ent rotation found in some elephant trunks ( Gahm et al .120061) . This 
may require an initially dynamic density field si nce such structures 
seem t o occur quite readily in the simulations of IGritschneder et aP 
(I2OIQI) . 

Our simulations (and those o f lHennev et al.l2009h suggest that 
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ionisation fronts with a strong perpendicular magnetic field should 
have clear observational consequences in both the morph ology of 
the re combination emission and of the dense gas. Henne v et alJ 
(I2OO9I) note that while the Orion Bar has a similar linear shape 
to the ribbon-like structures seen here, it is associated with the 
main ionisation/dissociation front itself rather than being an over- 
density within the H II region. The extended photoionisation front 
an d photon-dominated re gion seen edge-on in Ml 7 was modelled 
bv lPellegrini et al. j 20071) as bein^ suppor ted by magnetic pressure 
(observations bv' Brogan & Trolan3l2()0ll suggest a LOS magnetic 
field B > 100 /iG). Unfortunately much of this H II region is 
heavily obscured at optical wavelengths and so straightforward in- 
spection of Ha images is not particularl y revealing. In the exten- 
sive Hq survey of the Carina nebula by ISmith. Ballv & WalbornI 
(l2010ah . only the elephant trunk in 'Pos. 23' of fig. 1 shows ev- 
idence of a bright ridge in front of the trunk, suggesting that 
such features are not com mon. A clearer example is seen in NGC 
6357 teohigas et al.ll2004l) where a prominent bright ridge in front 
of an elephant trunk is the brightest Ha structure in the field of 
view. A similar structure is found in front of a massive pillar in 
NGC 3603 (Brandner etal. 2000). In both of these observations, 
however, the adjacent massive star cluster is expected to drive out- 
flowing gas and so it is unclear if the bright ridges are ram-pressure 
or magnetically confined. 

The larger simulation domains used here (compared to ML 10) 
allow us to follow the evolution of the dense gas for significantly 
longer, up to 700 kyr in some simulations. We find that in the ab- 
sence of dense gas further from the star, a pillar-like structure will 
be flattened into a cometary globule with a dense head and low 
density tail. To study this in more detail, clump configuration 2 
was set up with an extra clump a further parsec from the radiation 
source. It was found that the pillar-like stage of the evolution lasted 
^ 600 kyr in this model, until all of the dense gas reached the po- 
sition of the furthest clump. Subsequently a cometary morphology 
developed again. This experiment shows that while the pillar's life- 
time can be extended somewhat, it cannot survive indefinitely un- 
less there is a long dense filament pointing away from the radiation 
source. This supports the general picture of elephant trunks shown 
in fig. 15 of Smith et al. ( 2010b) as structures which form and dis- 
perse over about 10^ yr. With our non-gravitating simulations we 
cannot model the star formation w ith al so occurs, but the s imula- 
tions by ICritschneder et all feOlOh and iBisbas et all feOlOh show 
that the compression induced by RDI produces gravitationally un- 
stable fragments which would likely form stars. While it may 
seem unlikely that the three pillars in Ml 6 should all have origi- 
nated from elongated overdensities pointing back towards the ion- 
ising stars, recent o bservations from the Herschel observatory (e.g. 
iMolinari et al.ll2010 ) have shown that molecular clouds appear to 
have a distinctly filamentary struct ure. Simulations of MHD turbu- 
lence generated by colliding flows iBaneriee et al.l2009 ) show sim- 
ilar structures, an d simulations of small-scale n on-ideal isothermal 
MHD turbulence JPownes & O'Sullivanll2009l) also show signifi- 
cantly more linear structure compared to purely hydrodynamical 
models. 



7 CONCLUSIONS 

We have performed a series of R-MHD simulations of the pho- 
toionisation of dense clumps of gas and their evolution from pillar- 
like to cometary globule-like structures. Our results for the emis- 
sivity of ionised gas agree very well with those of Henney et al. 
([2009>) . showing that a dense, ionised, bar- shaped region standing 



off from the ionisation front is a generic feature of strongly magne- 
tised photoionisation in a clumpy medium for a perpendicular field 
orientation. This ridge can be as bright as, or even brighter than, 
the photoionisation front when observed in recombination radia- 
tion (e.g. Ha) and its presence or absence can be used as a diag- 
nostic of the strength of any large scale magnetic field which may 
be present. Bright ridges or ribbons are observed in som e H II re- 
gions (e.g. Brandner et al. 2000; Bohigas et al. 2004; S mith et al.l 
l2010ah . although they are not common. An overdense ridge could 
also be produced by ram-pressure confinement, and more detailed 
modelling is required to find observational signatures which could 
distinguish these different confinement mechan isms. 

Comparing to observations of Ml 6 (Hester et al ][l99i) there 
is no such ribbon or ridge, suggesting the ambient field measured 
by ISugitani et al.l 12007) is not dominant in the ionised gas. This 
conclusion is strengthened whe n we consider the magnetic field 
orientation observed in Ml 6 (Sugitani et "Zl bOOl fig. 9). The re- 
sults presented here show that both RDI and acceleration of clumps 
by the rocket effect tend to align the magnetic field in dense neu- 
tral gas with the radiation propagation direction. In our models 
a field configuration similar to the observed one is clearly seen 
when the initial field strength is |B| 20 /iG; the simulation with 
|B| ~ 50 /iG is consistent with observations, and the simulation 
with |B| ~ 160 /iG is not consistent. Our simulations thus suggest 
an ambient field strength of |B| < 50 /iG around the Ml 6 pillars. 

The morphology of the structures which develop due to RDI 
and the rocket effect is also affected by a strong magnetic field, 
partly due to shielding by the dense ionised ridge and partly by 
the effect of the field within the pillar or globule. Inspection of 
Ha images of elepha nt trunks and globules in the literature (e.g. 
Hester etal. 1996; S mith et al.L 2010a) suggests that the features 
seen in the strong field simulations are not common, although the 
uniform initial field configurations considered here are certainly 
somewhat artificial. Additionally many H II regions have signifi- 
cantly higher gas pressure than that modelled in our simulations, 
in which case the magnetic field must also be correspondingly 
stronger to dominate the dynamics. 

Comparing these results with our earlier simulations in ML 10, 
the larger simulation domains used here show that the pillar-like 
structures which form will ultimately evolve to cometary structures 
in the absence of dense gas further from the star. The lifetimes of 
pillars in our models are t < 500 kyr, although this depends sig- 
nificantly on the initial mass and concentration (and presumably 
velocity, cf. Gritschneder et al. 2010) of the dense gas clumps. 

F inally we emphasise, in agreement with previous au - 
thors t lWilliamsll2007l : lKrumholz et al.ll2007l : lHennev et al.ll2009l) . 
that a strong magnetic field has a very significant influence on the 
dynamics of the photoionisation process, and many of these effects 
should be easily observable. Given the difficulty of measuring the 
full 3D magnetic field in the ISM, comparison to detailed numeri- 
cal simulations such as these offers an indirect means to constrain 
the field strength and orientation in and around H II regions. 
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Figure Al. 2D shock tube results for the Brio & Wu (1988) adiabatic MHD test problem. The solid line shows the results using 10,000 grid cells on a ID grid, 
and the points show results in 2D on a domain with 200 x 200 cells at time t = 0.12. The panels show (from left) gas density, pressure, normal velocity and 
tangential field. Note that because all points in the 2D domain are plotted, the number of points apparently contained within a discontinuity is larger than the 
actual number of points which resolve it. 




Figure A2. Circularly Polarised Alfven Wave test. The left-hand plot shows the transverse magnetic field after advecting five times across the domain (initial 
amplitude 0.1) for models with resolutions A^^ = [16, 32, 64, 128]. The right-hand plot shows the LI error (see equation IaB as a function of resolution for 
this test, obtained by comparing the initial to the final solution at each resolution. Second order convergence is clearly demonstrated. 



APPENDIX A: TESTS OF ADIABATIC MHD 
ALGORITHMS 

Results from the following MHD test problems can be found at 
http : / /www .astro . uni-bonn . de/~ jmackey/ jmac/ 
together with the HD test problems. These adiabatic calcula- 
ti ons are run using dimensionless units. The shock-tube tests 
of iBrio & wJ (Il988h and lFalle et all (Il998h were performed in ID 
and on a 2D grid with the shock-normal at an ang 
to the grid axes, ensuring that the probl em is genuinely m ulti- 
dimensional. 2D Results are shown for the lBrio & Wul (Il988h test 
in Fig . lAll and are comparable to those obtained by Ipalle et aP 
(Il998h . The propagation of a polarised Alfven wa ve w a s also 
calculated in 2D, using the initi al con ditions from TothI 
and set up on a grid as in Stone et al ] (120081 Running this test 
at different resolutions shows the degree of numerical diffusion 
in the algorithm and also the rate of convergence of the solution. 
Results are shown in Fig. IA2[ comparison of the first panel to 
results from Stone et al.l (12008 ) shows that the third order ATHENA 
algorithm is significantly more accurate for this test at a given 
resolution, as expected. The second plot in this figure shows the 
LI error: 



for N cells, where 0^ is the reference state for each cell and is 
the approximate numerical solution. This verifies that the code is 
indeed second order accurate for continuously differentiable den- 
sity fields. 

^Results for the Orszag-Tang vortex problem ('Orsz ag & Tand 

[1971) were also very similar to those from the ATHENA code, with 
slightly more diffusion for smooth waves because of the lower or- 
der of accuracy. 

The advection of a magnetic field loop fe.g. lStone et al.|[2008b 
across a periodic domain is a challenging test for grid-based codes. 
When run with a non-zero velocity in the third dimension (vz 7^ 0) 
this test shows a weakness of the mixed-GLM divergence cleaning 
method - small divergence errors lead to small force errors and 
some of the initial magnetic field leaks unphysically into the z- 
direction 5 per cent). Results are shown in the first three panels 
of Fig. lA3l for the evolution of the magnetic pressure. The left panel 
shows the initial conditions; the centre panel shows the state at t = 
2 for a model with no advection (to show the numerical diffusion), 
and the right panel shows the state at t = 2 for a model advected 
twice across the domain. The time-decay of magnetic pressure for 
the advected models is shown in the right-most panel of Fig. I A3 1 
The constrained transport method employed in the Athena code, 
together with its third order accuracy, leads to less diffusion and a 
more sy mmetric solution to this problem than is obtained with our 
code fcf. lStone"et al...2008.) . but our results are comparable to those 
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Figure A3. Field Loop Advection test: the initial (1st panel, from left) and final magnetic pressure is shown, with static (2nd panel) and advected (3rd panel) 
models. Only the central half of the domain is shown. 15 contours are shown on a linear scale from the minimum to maximum values (dimensionless units). 
Some extra diffusion and distortion is apparent in the advected mo del. The d ecay of magnetic pressure (normalised to the initial value) over time is shown 
in the 4th panel for a model run using the van Albada slope-limiter Ivan Albada et al. 1982) . The three curves show the (small) effect of increasing artificial 
diffusion with the viscosity parameter iFalle et al .11 199 8) indicated. 



Time = 0.200 




APPENDIX B: BOUNDARY CONDITION FOR 
MIXED-GLM DIVERGENCE CLEANING ALGORITHM 



X X 

Figure A4. Contour plots of sHces through the gas density (left) and mag- 
netic pressure (right) at time t = 0.2 for the MHD blast wave test problem 
with /3o = 0.002. The slice is through the mid-plane z = 0.5. There are 30 
contours on a linear scale from the indicated minimum to maximum values. 
B ecause of the higher wave speeds here compared to the /3o = 0.2 model 
in I Stone & GardineJ i2009h . the fast magnetosonic waves seen in their re- 
sults have long since left the domain in this simulation. 



of other 2nd order astrophysics MHD code s. Tests also showed tha t 
for this test the van Albada slope-limiter I van Albada etaDll982h 
produced a noticeably superior solution than the commonly used 
MinMod limiter, probably because the latter is significantly more 
diffusive. 

The expansion of a multi-dimensional adiabatic blast wave 
provides a good test of a code's robustness in modelling 3D 
shocks and rarefactions. W e ran the same 3D problem as described 
in lStone & GardineJ (l2009l . §6.5), with the gas pressure in a sphere 
at centre of the domain set to 100 x the ambient pressure, and with 
a uniform magnetic field at 45° to the grid axes. Very similar re- 
sults were found for the case where the initial plasma beta parame- 
ter/3 = 0.2. Because of the robustness of the Dedner et al. ( 2002) 
algorithm, and presumably also because the integration algorithm 
used is slightly more diffusive than the ATHENA algorithm, our 
code can also simulate this test problem with a field 10 x stronger 
(/3 = 0.002), shown in Fig. lA4l This extra robustness is important 
for the strong field simulations run here. 



iDedner et al.l (l2002h suggest that for zero-gradient boundary con- 
ditions (BCs), using zero-gradient for the unphysical scalar field 
ip is sufficient. We have found, however, that more care is needed 
for magnetically dominated regions (plasma parameter /3 = 
Sirpg/B^ <^ 1) when gas is flowing subsonically near the bound- 
ary. In what follows x is the normal vector to the boundary, the 
computational domain is 'left' of the boundary, and off the domain 
is 'right' . Note that is not a physical variable so its gradient is 
not constrained. In fact if the field is constant across the bound- 
ary then the formal solution to the ID Riemann problem is that the 
flux, F[Bx] = (as it is for all ID Riemann problems). So ide- 
ally a boundary c ondition for t/ ^ shoul d be selected which ensures 
F[Ba:] = 0. From lDedner et alJ (|2002, eqs. 41,42) it is easy to cal- 
culate this: the zero-gradient condition states that — = so 
if additionally ip^ — —ip^ then F[Bx] = is obtained. 



A sequence of figures from a 2D slab- symmetric test simula- 
tion (this time using c.g.s. units) are shown in Fig. IB II to demon- 
strate the effectiveness of the new boundary condition. The simula- 
tion domains used have sizes [0.2, 0.1] pc (large) and [0.1, 0.1] pc 
(small) and the radiation source is at infinity in the — x direc- 
tion. A group of randomly located dense clumps are placed be- 
tween X = 0.02 pc and x = 0.08 pc on a uniform background 
density of nu = 100 cm~^. The initial magnetic field is B = 
[50, 1, 0] X ^/Att /J.G, and the initially constant gas pressure is set 
to give 13 = 0.017 (T = 1500K in the lowest density gas). Zero 
gradient BCs are enforced in the =bx directions, and periodic in 



Fig. IB II shows the gas pressure at time t = 12 kyr using cool- 
ing model C2. The problem with the zero-gradient BC for ip is 
very apparent, as is the dramatic improvement on switching to the 
= -^^ BC. The black regions are hiding cells which obtained 
negative pressures and were set to an artificial pressure floor value 
in order to allow the simulation to continue. Regarding the other 
boundaries, the x = 0.2 pc boundary only has waves moving per- 
pendicular to it and so there is negligible net flow across the bound- 
ary. The X = pc boundary has strong supersonic outflows, so the 
details of the BC for have little effect on it. 
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Figure Bl. Plots of gas pressure on a log scale from 10~^^- 
10~^ dyne .cm~^, at time t = 12kyr. Ionised gas has typical pressures 
of 10~^ dyne .cm"-^ and neutral gas 10"-*^^ dyne .cm~^. The top panel 
shows the full domain with = as the BC. The second panel shows 
the same, but with a simulation of only the left half of the domain over- 
plotted to highUght the boundary effects, again with = ipK The third 
panel is as the second, but using = —ip^ as the BC for both simulations, 
and the bottom panel shows the full domain with the = —ip^ BC. The 
simulations were run using cooling model C2. 



Name 




3^max 


2/min/ Zmin 


2/max/2^max 


cells 


Small 


2.0 


4.0 


-0.75 


0.75 


128 X 96^ 


Medium 


1.5 


4.5 


-1.125 


1.125 


192 X 1442 


Large 


1.5 


5.5 


-1.5 


1.5 


256 X 1922 


X-Long 


1.5 


6.5 


-1.5 


1.5 


320 X 192^ 


X-Wide 


1.5 


5.5 


-2.0 


2.0 


256^ 



Table CI. List of simulations run to test the boundary effects on the results 
of 3D MHD simulations of photoionisation. All tests were run with a 15 x 
a/Itt ~ 53 /j,G magnetic field in the y-z plane. Boundary positions are 
measured in parsecs in coordinates where the source is at [0, 0, 0], as in 
Tablen 
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Figure CI. (1) Neutral mass, (2) Position of leading edge of ionisation 
front, and (3) mass-weighted mean neutral gas x-velocity. Note that the 
large, extra-long and extra-wide simulations have almost converged for all 
plots. 



APPENDIX C: BOUNDARY EFFECTS AS A FUNCTION 
OF DOMAIN SIZE 

To study boundary effects on the results presented here, test simu- 
lations using clump configuration 1 were performed with increas- 
ingly larger domains, both in 2D and 3D with a medium perpendic- 
ular magnetic field (model R5 in Table |2]). The same physical reso- 
lution was used for each simulation; boundary positions for five 3D 
simulations are listed in Table ICll Various global quantities were 
measured within the volume common to all simulations. Neutral 
mass, mean neutral gas velocity, and position of the first neutral 
cell on the domain are plotted as a function of time in Fig. ICll 
The small and medium simulations are clearly significantly affected 
by the boundaries whereas the large model is essentially identical 
to the X-long and X-wide models, indicating that the results have 
converged at least in this limited subdomain of the simulation. The 
position of the leading edge of the ionisation front is unaffected by 
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the position of the boundaries, but the neutral gas mass and veloc- 
ity are significantly affected. It was found that flow of gas through 
the distant boundary was impeded in these R-MHD simulations, 
something which did not happen in R-HD simulations. 
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